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Abstract 

We introduce stochastic Discrete Laplacian Growth and consider its 
deterministic continuous version. These are reminiscent respectively 
to well-known Diffusion Limited Aggregation and Hele-Shaw free 
boundary problem for the interface propagation. We study correla- 
tion between stability of deterministic free-boundary problem and 
macroscopic fractal growth in the corresponding discrete problem. 
It turns out that fractal growth in the discrete problem is not influ- 
enced by stability of its deterministic version. Using this fact one 
can easily provide a qualitative analytic description of the Discrete 
Laplacian Growth. 

1 Introduction 

The Laplacian Growth or zero surface tension Hele-Shaw free-boundary problem 
( see e.g. pQ, [6]) describes time evolution of a domain f2 = f2(i) in the complex 
plane with z = x + iy, when the boundary <9f2 of the domain is driven by the 
gradient of a scalar field = cf>(x, y, t) (see Fig. [1]). The field is harmonic, except 
for several singular points ("sources" and "sinks"), in the exterior or the inte- 
rior of f2. The corresponding Hele-Shaw problems are referred as exterior and 
interior respectively. The field cj> vanishes at the domain boundary (" interface" ) 
4>(dfl) = 0. Well possedness and stability of the problem depends strongly on 
types of singularities. 

For instance, the exterior "sink"-driven Hele-Shaw problem is a problem 
of expanding a simply-connected domain whose boundary is driven by the 
field harmonic in exterior of fl. It is ill-posed and linearly unstable almost for 
any initial conditions. On the contrary, the exterior problem for contracting Q 
(which is a "source" -driven time reversal of the above problem) is linearly stable 
and well posed. 

A discrete, stochastic version of the exterior Hele-Shaw problem (often called 
Diffusion Limited Aggregation or DLA, [2], [3], [7]) describes formation of a clus- 
ter of particles on two dimensional (e.g. square) lattice. No more than one 
particle can occupy a lattice cite. The particles stick one by one to the cluster. 
The probability for a particle to occupy a given (next to the cluster) cite is 
proportional to the value of a lattice harmonic field at that cite. The field has 
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logarithmic "sink" -type singularity at infinity. It vanishes on the cluster and is 
updated after each cluster increment. 

The cluster in such a discrete problem is a fractal and one may think that 
instability and ill-possedness of the deterministic continuous version of such a 
discrete model (i.e. exterior Hele-Shaw free-boundary flow for expanding domain 
fi) is a manifestation of this fact. In other words, the fractal interface of the 
discrete problem may not be approximated by an analytic boundary curve of 
its deterministic counterpart. 

In this article we consider a "mix" of interior and exterior problems for 
simply-connected bounded domain, when the interface is driven by the scalar 
field 4> which is harmonic almost everywhere except for the interface <90 itself, 
where <f> vanishes, and two logarithmic singular points. One of these points is 
placed at z — and other at z — oo. In such systems the growth depends on 
a parameter A which is a measure of the "mix" between the exterior unstable 
(A = 0) and the interior stable (A = 1) processes, with A = 1/2 being a "neutral 
stability" growth. 

One may expect that the fractal formations are not visible on macroscopic 
scale for discrete version of such models when A > 1/2 and when the interface 
is linearly stable. This turns out to be true only for extreme stability point 
A = 1 , while the macroscopic fractal formations are present for any A < 1 in the 
discrete model. 

Note that the Discrete Laplacian Growth (DLG) introduced here differs from 
the Diffusion Limited Aggregation (DLA) by simultaneous consideration of both 
exterior and interior processes on the lattice defined in a similar way. In the case 
of DLA one considers a lattice cluster and a complement to it, while in DLG the 
lattice is divided into interior domain, discrete boundary and exterior domain. 
Although the law of the cluster growth in pure exterior limit (A = 0) of DLG 
is locally different from that of DLA, both models have the same continuous 
version and belong essentially to the same class. 

In the next section we describe, in details, the continuous version of the 
growth processes under consideration and perform its linear stability analysis. 

2 Continuous Model 

Let, for simplicity, f2 be a simply connected, bounded domain in the z = x + iy 
plane with the point z — inside the domain. We denote by <j) + the field 
(f>(x,y,t) in the interior and by in the exterior of the domain respectively 
(see Fig. [T]). Field (j>± is harmonic in the interior/exterior of f2 except for points 
z = and z — oo 



Acf)± = 0, z £ {0,oo, dfl} . 
It is continuous across the boundary and vanishes on it 



(1) 



<t>±{z £dfl) = 



(2) 
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C,(a,t+dt) 



'(b,t+dt) 



Figure 1: Problem setting (left) and the domain area increment along the bound- 
ary segment a < I < b during the time interval (t, t + dt) (right) 



The field has logarithmic singularities at z — and z = oo 

<j> + -> ^log|z|, z^O, 
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log \z\, z — > oo 



where A± are constants. 

Consider the situation when the boundary dynamics is governed by gradients 
of <f>± 

d<j> + d<j>- 

v n = a + — h a_— — , 

On On 

where a± are constants, v n denotes the normal velocity of the boundary 

v n = Re(ndz/dt), z — x — iy, z £ 90, 

n = n x + in y , \n\ = 1 stands for the exterior normal to <90 and d/dn denotes 
the normal derivative at the boundary. 

Rescaling the harmonic fields <j>± and the time variable t 



t 



A 



a+A. 



-t, 



where 



A 



A + a + 



A- a- + A+a+ ' 
we rewrite the above dynamic law for the boundary velocity as 

d<f>- d(j)+ 
On on 



(3) 
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and field asymptotic now depend on a "stability" parameter A as 

-> log |«|, z^O, (j>- -> — log|z|, Z^OO. (4) 

Z7T Z7T 

Equations (|H2l3l4p together with the initial condition f2(t = 0) constitute a 
free-boundary, initial value problem for f2(i). 

Note, that the case A = corresponds to the exterior Hele-Shaw problem 
(with (j> + = 0) while A = 1 refers to the interior Hele-Shaw problem (with 
0_ = 0, respectively). 

We are interested in discrete stochastic processes of lattice cluster growth, 
corresponding to the deterministic problem (jTJJJ) ■ 

Consider first the case of cluster growing from a single particle at origin that 
corresponds to the circular solution of deterministic model. It is easy to see 
that such solution is an expanding circle of radius r(t) = \ftfn and <j>± — <p± , 
where 

^+ } = ^(log kl - ^ log -), |z|<r(t), 

Z7T Z 7T 

(o) = l^A (logN _l log * )) N>r(t) 

Z7T Z 7T 

Let us now study linear stability of this solution, considering small perturbations 
of the circle \z\ = r(t) by the fc-periodic function 

|js(0,t)| = r(t)+eA k} {t)sm(k9) 

where 9 = (0,27r) is a polar angle on the z-plane (z = \z\exp(i9)) and k is a 
positive integer. The fields <f>± satisfying conditions (JTJ) , (JD) are of the form 

4>± = 4>± ] + £0 (fe) {t)\z/r(t)\ ±k sm{k9). 

Substituting this and previous equations to ^ and (O in the first order of 
e-series we get 

dr[k) 1-2A, (k) 

— — = — kr w (5) 

dt 2t w 

Since k > 0, the perturbations grow in time when A < 1/2 and vanish in time 
when A > 1/2, with A = 1/2 being the "neutral stability" point. 

In this article we consider continuous systems labelled by stability parameter 
A in the range < A < 1. The A = case (exterior Hele-Shaw problem for 
expanding interior domain) is unstable and ill-posed problem, whose discrete 
version manifests fractal growth. On the other extreme of the interval, at A = 1 , 
we have stable and well-posed interior expansion problem. 

Now it is natural to ask the question: Whether the stability of interior con- 
tinuous problems (i.e. those at A > 1/2) may depress the macroscopic (i.e. 
visible in continuous limit) fractal growth of boundaries of corresponding dis- 
crete models? 

To study this question, we are to introduce discrete analogue of the above 
problem ([HI]). 
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3 Discrete Model 



One may think of the free boundary problem (jTj4j) as that of dynamics of an 
ideal conducting contour dQ in two-dimensional electric field. The field is cre- 
ated by a Coulomb charge of value — A at z — and unit charge distributed 
with linear density cr(l,t) along the contour. The contour is "an ideal conduc- 
tor" , and this means that the potential (at fixed t) is constant along dil. Here 
/ stands for natural parameter along the contour {dfl : z = ((l,t), \d£(l)\ = dl} 
and § dn a(l,t)dl = 1. The harmonic field <fi(x,y,t) is (modulo coordinate inde- 
pendent function of t) the electrostatic potential created by the above charges 



^ log \ z \ + vQ, *) i°g I* - CG, t)\di 



The gradient of the potential is the electric vector field, that has jump of mag- 
nitude a (I) across the boundary. Therefore from it follows that the normal 
velocity of the contour equals its linear density 

v„(l,t) = a(l,t) 

This and the previous equations reduce the free-boundary problem on the plane 
to a dynamical problem on the contour only. It is easy to see that the density 
a is non-negative, and the interior domain fl expands in time. 

The domain area increment dP( a m along the interface segment a < I < b 
during the time dt (see Fig. Q]) 

dP(g, b) 

dt 



v n (l)dl= I a{l)dl = q {afi) (6) 



is proportional to the electric charge q( a ,b) concentrated on this segment. This 
fact suggests to consider the boundary charge as the cluster increment proba- 
bility in discrete analogue of the above model. 

Let us now take a square lattice, paint elementary square cells in white 
color and label centers of squares by a pair of integers (m, n) that are discrete 
coordinates along the x and y directions. Consider the following stochastic 
process: At the first step, the (0,0) square is painted in black. At the next step 
we paint in black one of four cites that are next-neighbors of this smallest cluster 
e.t.c. (see Figure[5|). This procedure is continued by coloring, at each step, any 
randomly chosen white next-neighbors of the cluster (i.e. adding randomly a 
"particle" to the cluster) with the probability described below. 

We use the same notation ft for the cluster as we used for the continuous 
domain. The cluster boundary, which consists of all white next-neighbors of the 
cluster, is denoted, by analogy with continuous case, by dfl. 

At each step of the process we can define a lattice function m ,„ which 
satisfies the difference equation 

<f>n-l,m + 4>n+\,m + 4>n,m-l + K,m+1 ~ #n,m = -A5„ )0 ^m,0, {n, m) g dfl (7) 
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Figure 2: An example of the cluster growth. Black squares ("particles") belong 
to the cluster 0, white squares with gray sides form the cluster boundary d£l. 
Black dots denote charged (boundary) cites and gray dots are uncharged lattice 
cites. 



everywhere except for the boundary, where it vanishes 

<t>n,m = 0, (n,ro)€0fi. (8) 
Its asymptotics is like follows 

— > 1 , log(m 2 + n 2 ), m 2 + n 2 — > oo (9) 

47T 

The left hand side of is a lattice Laplace operator, and 6 on the right hand 
side stands for Kronecker delta symbol. Equations (JTHH]) are lattice analogues 

of dull!. 

As in the continuous case, the potential <j) n . m can be expressed in terms of 
charges, that are now placed at the boundary cites (see Fig. [2]). 

Introducing the lattice Coulomb potential (or the Green function) G n ^ rn , a 
such that 

- 4G n , m = 5 n fiS m ,o, — °o < n, m < oo 

and 

G n m — > — log(m 2 + n 2 ), m 2 + n 2 — » oo 

47T 

one expresses n , m in terms of q,;^, (i, j) S 9il as 

where C is a constant. 

It's easy to see that the boundary charges are nonnegative % j > 0. From 
(J7J), © and by the definition of the Green function it follows that 

E '/<..- ( 10 ) 
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By direct analogy with continuous case (c.f. ©) we now consider as a 
probability of the square of the boundary dtt to join the cluster fl. 

Therefore, at each step of the growth process we have to solve the linear 
algebraic system consisting of equation (JTUJ) and 

'!. ,■(■, .-j y \C,j I C 0. (i,j)€dtt (11) 

for unknowns qij, € dtt and C, then repaint one of the boundary squares 
in the black color with the probability q^j. 

Note, that the lattice function for boundary cites satisfies the following equa- 
tion 

4>i-l,j + <Pi+lj + <i>j-l,i + 4>j+l,i = ftj> e dfl 

and the cluster increment rule can be also reformulated in terms of the potential 
4> as follows: Probability of the boundary cite to join the cluster equals to 
the sum of potentials on the next-neighbors of that cite. 

4 Fractal and Continuous Growth 

Since in the case of pure exterior A = problem the model under consideration 
differs from other discretizations of the Laplacian Growth by local details only, 
we have to expect the qualitative behavior to be similar to that of Diffusion 
Limited Aggregation. Indeed, for A — one observes the pure fractal growth 
(see Figure ■ Numerical calculations give the following estimate for Hausdorff 
dimension of fractal boundary of a tree-like cluster (see Appendix for details of 
numerical calculations) 

d= 1.64 ±0.02 (12) 

Note that in the A = case the dimension of the boundary coincides with that 
of the cluster. 

In the case of the pure interior problem A = 1, the Discrete Laplacian Growth 
shows dynamics, close to that of its deterministic continuous counterpart: The 
cluster boundary is not a fractal, and tends to solution of the continuous problem 
as its size increases (see Figure[3]). Such a behavior was expected due to stability 
of the interior continuous problem for expanding cluster. 

Since the continuous problem (HE]) is stable for A > 1/2 (c.f. ©) one would 
expect similar (close to deterministic) behavior of Discrete Laplacian Growth 
for such A. This turns out not to be so: Instead, one observes separation of the 
growth in two fractal components, one of which is quasi-regular and the other 
is macroscopic, i.e. visible in continuous limit, constituting a finite fraction of 
the cluster. In the range < A < 1, a big (grown from a single particle) cluster 
consists of a quasi-circular center and branches going out of it (see Figure 2]) . 
The Hausdorff dimension of such cluster boundary equals the same d as in (fT2^) . 

One can interpret such a behavior in the following way: The discrete growth 
can be viewed as a competition between the averaged process tending to a 
continuous limit, and probabilistic fluctuations. Such fluctuations drive the 
system beyond the linear stability region when A < 1. 
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Figure 3: Clusters (left) and their boundaries (right) for the pure exterior (A = 
1) problem (top) and the pure interior problem A = (bottom). To visualize 
the time dynamics, we prescribe colors to cluster particles. The color ranges 
from green to blue, depending on time (or step) at which the particle joined the 
cluster (green for the initial step, blue for the final step) 
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5 Qualitative Analysis 



One might try to explain the above fractal properties of growth in the stabil- 
ity region 1/2 < A < 1 by the fact that the cluster evolution considered so 
far starts rather with a single particle than with quasi-continuous macroscopic 
initial conditions. And indeed, to describe properly the continuous limit we 
have to start with a cluster of dimension 2, containing N — > oo particles. It is 
also necessary that, after rescaling the lattice spacing (or cluster linear size) by 
the factor iV -1 / 2 , the initial boundary dfl(t — 0) tends to an analytic curve as 
N — > oo. The time, which is the step number in the discrete model, is rescaled 
by the factor 1/N. 

Such growth process remains, nevertheless, a superposition of smooth (an- 
alytic) and macroscopic fractal components in the limit N — > oo: The quasi- 
analytic part of the boundary dCl(t) evolves as in the pure interior problem, 
while the quasi-analytic growth is superposed with the fractal growth of the 
pure exterior problem. In other words, the "mixed" interior-exterior problem 
separates into interior and exterior parts in zero-lattice spacing limit. 

It is easy to find the rate of growth of quasi-analytic to fractal components 
provided the separation always takes place for < A < 1: The tree- like fractal 
part mainly grows on the tips of the "branches", since the exterior field </>_ 
is screened out in fjords between the branches, while the interior field is 
screened out inside the branches. Therefore, in zero-lattice spacing limit, the 
interior (quasi-analytic) part of the boundary evolves like there were only one 
logarithmic singularity 

<j>+ ^ -^\og\z\, z^O, and = 0, 

while the fractal part of the curve grows as if 

4>- — > — - — log \z\, z — > oo, and </>+ = 0. 

The ratio of growth between the fractal and quasi-analytic parts (in number of 
particles per unit of time) is 

(1-A)/A. 

It is also straightforward to show that the separation always takes place for 

< A < 1. Let's consider the case when A is close to 1, since the separation of the 
fractal component here implies the separation for a smaller positive A. Suppose 
again that initial cluster has a big size N — > oo and its rescaled boundary is 
close to an analytic curve, so the number of the boundary cites is of order \/N. 

Consider now a perturbation (fluctuation) of the boundary. Suppose that the 
fluctuation linear size (in lattice spacings) is I. For A close to 1, the charge of the 
particles on the fluctuation tip (i.e. most distant from the quasi-analytic part 
of the boundary cites of the fluctuation) is of order (1 — A)/ 1- - / 2 /V~N, where 

1 < D < 2. For instance, D = 1 for one-dimensional, "crack"-like fluctuations of 
length I, while D — 2 for a "bump" -like 2-dimensional fluctuations of diameter 
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I. Charge of particles on the quasi- analytic part of the boundary is of order 
1 / vN. A fluctuation tend to grow (rather than collapse) when the probability 
for a particle to join the cluster at the tip of fluctuation exceeds the probability 
for the quasi-analytic part of the boundary. Since the probability is proportional 
to the charge, from the above discussion it follows that fluctuations tend to 
grow when their linear size / exceeds some critical value which is proportional 

_ 2 

to (1 — A) 2 -° . This value is asymptotically independent of cluster size N. 

Since the critical size of fluctuation is TV-independent, the possible number of 
critical fluctuations along the boundary is proportional to the boundary length 
>/N. It then follows that starting from the close to analytic boundary, a big 
cluster develops critical fluctuations in the number of steps f(\)y/~N, where 
/(A) is some function of A only. This number of steps corresponds to the time 
t = f(X)/V~N in the scaling limit. 

Therefore, in the N — > oo scaling limit, the cluster develops critical fluctua- 
tions in zero time and the separation of growth into the fractal and continuous 
components takes place immediately if < A < 1. 

6 Conclusions 

In the present work we have introduced the Discrete Laplacian Growth (DLG) 
model and studied correlation between the stability of its continuous version and 
the fractal growth of boundary. The stability, turns out, does not guarantee an 
absence of the fractal growth. The growth process separates into superposition 
of "continuous" and fractal component. In other words, in the N — > oo scaling 
limit, the discrete growth does not tend to solution of continuous model even if 
the later is stable and well-posed. 

There remain a few questions to address: 

1) Since, in N — > oo limit, DLG separates in two independent processes, one 
may think about implications of such a separation in the continuous model. For 
instance, whether this separation can be traced (in some, possibly "hidden", 
form) in continuous model or it is just a consequence of the discretization? 

2) It is also known that both, exterior and interior problems possess an 
infinite number of conserved quantities (harmonic moments) in the scaling N — > 
oo limit [51 [5] . Does this imply (in view of the separation) existence of an infinite 
number of conserved quantities in "mixed" models? 

3) In this article we presented numerical estimate flT2"j) for the fractal dimen- 
sion of DLG. Although DLG differs from the Diffusion Limited Aggregation 
(DLA) only locally, their fractal dimensions do not coincide. This result is not 
unexpected, since dimension of DLA is sensitive even to the lattice type (see 
e.g. [5]) and can not be considered as a general characteristics of a class of 
models. It would be interesting to understand which quantities are universal, 
i.e. independent of local details of a discrete scheme. 
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7 Appendix 



Numerical simulations of the Discrete Laplacian Growth have been performed 
by updating the discrete boundary (presented by a set of pairs of integer coordi- 
nates («(1), j(l)),(i(2), j(2)), . . . , (i(K) , j (K)) , where K is the current boundary 
length) with the probability given by solutions of linear system (llOlllj) . In this 
method an explicit solution of Laplace equation at each step is not needed, since 
the coefficients of Eq. (jTTJ) are defined by the coordinates of the boundary cites 
only. 

The Green function (or lattice Coulomb potential) G n ^ rn can be computed 
in several ways. To speed up the calculations, we have used the following ap- 
proximation for the lattice Green function 

1 4 i A n 2 2 

1 2 2 m + n — om n 
G m , n = — log m + n ) — 2^3- 

when the point (n, m) is in the exterior of the square —L < n < L, —L < 
77i < L, and numerical solution for G„. m in the interior of this square (with 
the function value at the square boundary given by the above approximation). 
In this approach, the Green function is represented numerically as (computed 
once) 2L + 1 x 2L + 1 table extended by the above equation. This approximation 
is very precise: for example, for L = 100 the maximum deviation from the exact 
lattice Green function is of order 10~ 10 . 

We split solution of the system (jlOllip in two parts: 



K 

E 

k' = l k' = l 



i(k), (k) 



r , — n„ cx w n r< — \ \ S fc=i fffc 1 
q k - -Cq k + Aq k , C — w 



Z^fe=i <ik 

where q k stands for qi( k ),j{k)- The K x K matrix G i ( k )-i( k i),j(k)-j(k') is sym- 
metric and positive definite, so we made use of the rapidly converging conjugate 
gradient method [3] for the numerical solution. 
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